NASTRAN THERMAL ANALYZER 


A GENERAL PURPOSE FINITE-ELEMENT HEAT TRANSFER COMPUTER PROGRAM 

Hwa-Ping Lee and James B. Mason 
NASA Goddard Space Flight Center 


SUMMARY 

A general purpose heat transfer analysis capability based on the finite-element 
method has been added to the NASTRAN system. The program not only can render tem- 
perature distributions in solids subjected to various thermal boundary conditions, includ- 
ing effects of diffuse-gray thermal radiation, but is fully compatible in capacity and in 
the finite-element model representation with that of its structural counterpart in the 
NASTRAN system. The development history of the finite-element approach for determin- 
ing temperatures is summarized. The scope of analysis capability, program structure, 
features, and limitations are given with the objective of providing NASTRAN users with 
an overall view of the NASTRAN Thermal Analyzer . 

INTRODUCTION 

The NASTRAN Thermal Analyzer is a general purpose heat transfer computer pro- 
gram based on the finite-element approach. It is fully capable of rendering temperature 
distributions and heat flows in solids subjected to various boundary conditions, including 
the modes of convection and radiation, in both steady-state and transient problems. This 
computer program has been achieved by the addition of new capabilities to NASTRAN in- 
cluding new elements and solution algorithms. The resulting Thermal Analyzer is both 
self-contained and user-oriented but remains totally compatible in capacity and in the 
finite-element model representation with that of its structural counterpart in the NAS- 
TRAN system. Therefore, thermal and structural analyses of large and complex struc- 
tural configurations utilizing a unified finite-element representation have become a 
reality. 

The purpose of this paper is to provide the NASTRAN users with an overall view of 
the NASTRAN Thermal Analyzer describing its scope, program structure, features and 
limitations. The practical needs which motivated the development of this thermal analy- 
sis capability are given. A brief survey of early works with regard to heat transfer 
analysis by the finite- element method is included. Also described are studies conducted 
at GSFC which investigated element behaviors, obtainable accuracy and solution efficiency 
of various schemes used in conjunction with the finite-element method, the feasibility of 
adding thermal analysis capability to the NASTRAN system, and the various approaches 
for treating problems with emphasis on radiative exchanges . 
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BACKGROUND OF DEVELOPMENT 


The main objective of a research and development program at GSFC, titled the 
Structural- Thermal-Optical Program (STOP) (ref. 1), is the devlopment of analytical 
methods and procedures which can yield reliable predictions of thermal deformation and 
optical degradation in an orbiting spacecraft (ref. 2), or in a simulated environment. 

This multi-disciplinary project is intended to render analytical services to projects such 
as LST (Large Space Telescope), SAS (Small Astronomy Satellite), EOS (Earth Observa- 
tory Satellite), etc. To predict alignment or optical performance of a large spacecraft- 
telescope system, NASTRAN is relied upon for structural analysis to compute thermally 
induced deformations and stresses. For a reliable structural solution, the thermoelas- 
tically uncoupled structural analysis requires accurate temperature data as an input to 
the NASTRAN structural model. Prior to the existence of the NASTRAN Thermal Analy- 
zer, available general purpose thermal analysis computer programs were designed on the 
basis of the lumped-node thermal balance method (e.g. ref. 3). They were not only -limi- 
ted in capacity but seriously handicapped by incompatibilities arising from the model 
representations inherent in the two distinct approaches. The inter model transfer of 
temperature data was found to necessitate extensive interpolation and extrapolation. 

This extra work proved not only a tedious and time-consuming process but also resulted 
in compromised solution accuracy. To minimize such an interface obstacle, the STOP 
project undertook the development of a general purpose finite- element heat transfer 
computer program. 

Preliminary studies were then conducted at GSFC to investigate the feasibility of 
achieving such a computer program as an integrated part of the NASTRAN system. 

When this task began, the theoretical foundation of the finite-element approach had been 
established for both steady-state and transient thermal field problems (ref. 4-6). Efforts 
aimed at broadening its scope were found in literature such as an extension to axis- 
symmetric problems (refs. 7, 8), a consideration of inhomogeneous material using higher 
order elements (ref. 9), a demonstration of procedure enabling the achievement of an 
efficient solution by reducing the order of the set of equations (ref. 10), and many special 
applications (ref. 11-13). All of these studies, however, were limited to thermal conduc- 
tion with linear boundary conditions. Since radiation is a dominant mode of the heat 
transfer process in space- oriented applications and introduces a fourth-power nonlinear 
temperature term in the boundary conditions, the in-house studies were directed princi- 
pally to include these nonlinear radiation effects. Other investigations into solution 
feasibility, accuracy and efficiency, and element behaviors in combined modes of heat 
transfer analysis were also undertaken. The in-house studies consisting of two parallel 
efforts were: 

(1) the structural version of NASTRAN was employed directly to solve heat trans- 
fer problems by utilizing mathematical analogies together with manipulations of 
available elements and solution routines (refs. 14, 15). The concept of treating 
the radiative fluxes as nonlinear loads was demonstrated successfully for simple 
radiative problems. The solution scheme, however, was restricted in that all 
nonlinear radiation problems had to be solved using the direct transient 
integration algorithm available in the program. 
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(2) the derivation of heat elements including radiative effects from the governing 
equation and boundary conditions approached via a variational principle, and 
studies centered on heat elements with associated different solution schemes 
(ref. 16). Two distinct approaches to nonlinear radiative problems were investi- 
gated and studies of radiative exchanges between element surfaces of known 
and unknown temperatures were included. A direct steady-state solution via the 
consistent linearization method and using an iterative process proved to be able 
to yield a solution more efficient than using a transient route. The other ap- 
proach, the direct energy distribution method, which had been illustrated inde- 
pendently in a simple problem (ref. 17), in fact, shares the same theoretical 
basis as that employed by a direct use of the structural version of NASTRAN 
(ref. 14, 15). A prototype computer program based on the direct energy distri- 
bution method was coded (ref. 18), and studies of element behaviors and solution 
accuracy were also conducted. 

While the NASTRAN System Management Office at Langley Research Center planned 
the extension of the NASTRAN system to include linear conductive heat transfer analysis 
capability only, GSFC-STOP was seeking the implementation of a full-fledged finite- 
element heat transfer computer program. The software capability was finally developed 
and implemented by the MacNeal-Schwandler Corporation under a subcontract from Bell 
Aerospace Company. It must be stressed, however, that a cooperative financial and 
technical effort between these two NASA centers made possible the emergence of this 
vital new capability. 


SCOPE OF ANALYSIS CAPABILITY 

The NASTRAN Thermal Analyzer is capable of solving both steady- state and trans- 
ient heat transfer problems in large size systems, of arbitrary configuration. The gov- 
erning heat equation together with the boundary conditions may be expressed, following 
the finite-element technique, in matrix form as 


[b] {d} + m (u> = (p> + m 


a) 


where 

{u} = vector of grid point temperatures 
{u} = rate change of temperature vector 
OK] = thermal conductance matrix 
[B] = thermal capacitance matrix 
{P} = vector of grid point heat fluxes 



{N} = vector of temperature- dependent nonlinear heat fluxes 

With appropriate interpretation of the matrices appearing in Eq. (1), three basic types of 
problem can be identified. They are distinguished by the numerical solution algorithms 
required for their solutions. 

(1) Linear Steady-State Problems 

Letting (u) = (N) =(0) and M = CKjl , Eq. (1) is reduced to a linear steady-state 
equation of the form 


og {u> = {p> 

where CkJ is a constant thermal conductance matrix. This equation is analogous to the 
basic linear static analysis of the structural version and is treated by the solution 
algorithm present in Rigid Format 1, (Ref. 19). 

(2) Nonlinear Steady-State Problems 

Letting (u> = (o) and Ck3 =[K 2 ], Eq. (1) has the form 


CK 2 ] {u} = {P} + (N> 

where [K 2 1 may be temperature-dependent. The permitted nonlinearities may arise from 
temperature-dependent conduction and convection properties as well as diffuse radiation. 
A new iterative solution algorithm, which includes test for convergence, has been added 
to the program for this type of problem. 

(3) Linear and Nonlinear Transient Problems 

Equation (1) directly represents the general transient heat transfer problem in 
which {P} is allowed to be a time-dependent heat flux vector. The permitted nonlineari- 
ties may arise from the effects of radiation or from user specified nonlinear elements. 
The new solution algorithm accommodated for this type of problem is implicit and is a 
combination of forward and backward differencing embracing a free parameter that 
allows the user to select the amount of each. While a default value exists for the param- 
eter which assures stability in linear problems, an override option is provided which al- 
lows the user to trade stability for efficiency. An option is also available which permits 
the user to linearize the effects of radiation. 

PROGRAM STRUCTURE 

The NASTRAN Thermal Analyzer has been designed to perform thermal analysis 
utilizing input and output formats compatible with those of the structural version. This 
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capability has been achieved by making maximum use of existing elements and system 
capabilities needed to satisfy the unique requirements posed by thermal problems. All 
input quantities as well as output displays are physically meaningful to users in the 
thermal field. 

Features such as nonlinear and scalar elements, multipoint constraints, transfer 
functions, DMAP, etc. are automatically inherited from NASTRAN. The user, therefore, 
is provided with a powerful tool to include effects not normally modeled by other ele- 
ments as long as those effects can be described by zero or first order system of equations 

The essential points which characterize this thermal analyzer are summarized as 
follows: 

(1) Geometry Description 

The body to be analyzed is idealized as an assemblage of appropriate finite 
elements. NASTRAN grid, scalar, and extra points remain valid for use in this descrip- 
tion. However, only one degree-of- freedom is associated with a grid point because of 
the nature of the scalar temperature field problem. 

(2) Types of Elements 

The program contains elements in three general categories: 

(a) Heat conduction elements — The constant gradient line, triangle and tetra- 
hedra are the basic elements and utilize linear temperature gradients in 
one, two and three dimensions, respectively. All elements share common 
descriptions with their structural counterparts and are summarized in 
Table 1. Quadrilateral elements are composed of overlapping triangles, 
and wedges and hexahedra from sub-tetrahedra. Solid-of-revolution ele- 
ments of triangular and trapezoidal cross-section are available for analyz- 
ing axisymmetric problems. Lumped capacitance matrices are used with 
all conduction elements to account for the effect of heat storage. 

(b) Boundary elements — XHBDY elements are used to define surface heat flux 
input and to describe convective and radiative exchanges between boundary 
elements. More than one boundary element can be connected to any one 
heat conduction element to account for multiple boundary heat inputs, in- 
cluding convection and radiation. A special cylindrical element with an 
arbitrary user specified elliptic cross-section is included for convenience 
in treating vector heat inputs. With this one dimensional element the ef- 
fects of directional radiant heat inputs are automatically included in that 
peripheral energy inputs over the lateral surfaces are integrated internally. 

(c) Special elements - Several element types are available for the indirect in- 
clusion of effects which cannot be otherwise modeled if only heat conduction 



elements together with boundary elements were used. These special ele- 
ments may be employed to enhance the analysis capability of the program. 
Scalar elements, nonlinear elements, transfer functions, and direct matrix 
inputs are included in this category. In heat transfer problems, scalar 
spring elements are analogous to thermal conductors and scalar dampers 
to thermal capacitance. 

(3) Material Properties 

Isotropic and anisotropic material behaviors are included. The treatment of 
nonlinearities arising from temperature-dependent thermal conductivity and convective 
film coefficient requires an iterative solution algorithm and this has been automated for 
steady-state problems only. The process involves the supply of temperature-dependent 
data in the form of tabulated functions which are interrogated at the beginning of each 
iterative step, hi transient cases, however, the user must rely on the use of nonlinear 
elements to treat the nonlinear effects as additional thermal loads that are evaulated at 
the previous integration time step. 

(4) Constraints and Partitions 

Constraint and partitioning features of the structural version remain valid. * 
Single-point constraints are used for the specification of prescribed temperatures, and 
multipoint constraints are used to describe known temperature dependency between tem- 
perature degrees of freedom. The omitted degree-of-freedom capability employs the 
well-known .Guyan reduction technique to reduce the problem size of the solution set. 

(5) Boundary Conditions, Initial Conditions, and Thermal Loading 

Convective exchange along the boundary can be specified between a surface and 
an ambient zone of known temperature, or between two or more surfaces or zones of 
unknown temperatures. Constant temperatures are specified directly by using constraints 
while temperatures which are arbitrary functions of time are specified indirectly by 
using simple modeling procedures which avoid unnecessary retriangularization when 
solving transient problems. Arbitrary initial temperature distributions can be specified 
in transient analysis. 

Several options are available to users for the specification of thermal flux in- 
put. Steady or time varying scalar heat flux can be described at the element and/or grid 
point level. Vector heat flux, such as that of radiant flux from a distant source, is 
described by specifying the flux intensity and vector direction. Both the flux intensity 
and vector direction may be time dependent and this allows, for example, the automated 
accommodation of rotating bodies in a vector heat flux field. Volumetric heat generation 
can be specified at the element level. 

Required input data for diffuse-gray thermal radiation problems consists of the 
Stefan- Boltzmann's constant, reference absolute temperature, thermo-physical properties 
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and an array of exchange coefficients together with a list which identifies the elements 
that are radiatively interacting. The array of radiative exchange coefficients AF is a sym- 
metric matrix according to the reciprocity rule. It is the product of the emitting surface 
area A and the geometric view factor F between the emitting and receiving surfaces. 

The Thermal Analyzer accepts radiation data via card or tape. 

(6) Graphics Capability 

The structural plotting feature for data checkout is available. Included options 
are orthotropic, perspective or stereoscopic projection capabilities. Time-history data 
of element heat flux, thermal loads, and temperatures at grid points can all be graphically 
shown in x-y plots. 

(7) Integrated Thermo- Structural Analysis 

The thermo-elastieally uncoupled thermal and structural analyses are performed 
in two passes through the NASTRAN System, but may be made to appear as a single con- 
tinuous run by the use of computer system control language. In the ease of transient 
analysis, temperature distributions computed by the NASTRAN Thermal Analyzer are 
recorded on magnetic tape or punched cards at predetermined time intervals and, sub- 
sequently, employed for static structural analysis. Back-to-back thermal and structural 
analyses require that the grid point locations must be identical in both models on a 
point-to-point basis. 


A VIEW FACTOR GENERATION COMPUTER PROGRAM 

In computing temperatures involving radiative interchanges between surfaces, the 
geometric view factors between any two radiatively active element surfaces are neces- 
sary to form the exchange coefficients as input to the Thermal Analyzer. In an in-house 
STOP project effort, GSFC has developed an IBM-360 program named ’'VIEW" which 
computes the view factors and the required exchange coefficients between radiating 
boundary elements. VIEW has compatible input-output formats with the Thermal 
Analyzer and possesses other programming features similar to those of NASTRAN. A 
detailed description of this program is presented in an associated paper titled "VIEW - 
A Modification of the RAVFAC View Factor Program for Use with NASTRAN Level 15." 

FURTHER STUDY IN EFFECT OF VIEW FACTOR COMPUTATIONS 

In studying accuracy, one of the factors that might influence the result of solution 
involving radiative exchanges between elements is the method that computes the amount 
of net radiant energy on the element level which is then evenly distributed to its vertices. 
This approximation involves the computation of view factors which are computed on the 
element surface to element surface basis associated with the element temperature which 
is an average of those temperatures at the vertices. An alternate approach is to form an 
isothermal area of the temperature at the grid point by dividing the origin elements into 
subelements and assembling the subelements from the adjoining elements to that grid 



point while subelements are formed by connecting the centroid to each midpoint on the 
side of an element. The difference of the view factor results from these two approaches 
is evident in view of a demonstration of the view factor Fa-b between two shaded areas 
as shown in Figure 1. 

A direct computation of F A B gives 0.23344, the result obtained from a summation 
of one-third of the view factors computed on the basis of element surface to element 
surface is 0.12471, and the result obtained from summing up the view factors on the 
basis of subelement to subelement following the rule of the view factor algebra is 0.233989. 
It is to be noted, however, that these computations took into consideration the view 
factor alone. A study combining the effects of temperatures and view factors which in- 
volves a modification of a prototype computer program to assess the quantitative influence 
to the solution accuracy is in progress and will be reported separately. 

ILLUSTRATIVE PROBLEM 

Since the delivery of the IBM-360 NASTRAN Thermal Analyzer to GSFC in June of 
1972, the system check-out phase has proceeded. A problem was designed to demon- 
strate program capabilities and features including: inter-element and inter-program 
(with regard to VIEW program) compatibilities, coordinate transformations, combined 
thermal modes of operation, vector heat flux input description, and the application of the 
different solution algorithms. The problem is that of a fin-like bent plate, whose dimen- 
sions, thermophysical properties and finite-element model are shown in Figure 2 where 
the input parameters with the subscript t refer to the tube and those without any refer to 
the plate. The underside of the plate is insulated thermally, and the upper face is exposed 
to a directional heat flux of S. The flowing fluid has a temperature T in = 1200 K at the 
entrance and a linear temperature drop of 60 K across the tube which has a wall 
thickness r = 1 cm and an outer radius r 0 = 20 cm. Determine temperature distribu- 
tions in the plate for the following cases: 

(1) The upper face dissipates heat to the surroundings of T ai r = 300 K with a con- 
vective film coefficient h air = 1.135 W/cm 2 -°K. 

(2) A radiative enclosure of a temperature Too replaces the convective environment 
in the preceding case, and radiative interchanges between all surface elements 
are taken into consideration. The solution is approached by the direct nonlinear 
steady- state solution algorithm. 

(3) The plate has a uniform temperature at T 0 = 300 K initially. The temperature 
response of the plate is determined by the nonlinear transient solution algorithm. 
The process starts with an activation of the flow into the tube which has identi- 
cal tube-fluid conditions as previously described. 

In view of the length that would be required to embed the complete meaningful com- 
puter print out of this problem, which consists of detailed card descriptions of the con- 
trol and input data decks and the output results, solutions are not included in this report 
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but have been prepared in a separately bound copy (ref. 20) which will be made avail- 
able to the audience of this colloquium as well as to any interested reader who will make 
a request to the authors. 
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TABLE 1 

HEAT CONDUCTION ELEMENTS 


DIMENSION 

TYPE 

ELEMENTS 

1-D 

Linear 

BAR, R(|>D, C0NR0D, TUBE 

2-D 

Membrane 

TRMEM, TRIA1, TRIA2, 
QDMEM, QUAD1, QUAD2 

Solid of Revolution 

TRIARG, TRAPRG 

3-D 

Solid 

TETRA, WEDGE, HEXA1, HEXA2 



PLANE-A 



Figure 1. Sub-elements joined at two common vertices to form two isothermal radiative surfaces. 




QVECT 

S = 0.14W/cm 2 


SP0INT 


T^, = 77.78 K 



P t = 2.707 
C t = 0.898 
h t = 0.8 


0.8 


k = 0.161 W/cm-K 
P = 7.849 g/em 3 
C = 0.492 J/g-K 

h - 1.135 x 10' 3 W/cm 2 -K 
ct = e = 0.9 


Figure 2. Finite-element model of a bent fin-tube configuration with boundary conditions. 



